function U = poisson_solver_2d(mesh, pde)
A = assemble_matrix_2d(mesh, pde);
b = assemble_vector_2d(mesh, pde);
[A, b] = process_boundary_2d(A, b, mesh, pde);
U = A\b;
end

%% assemble_matrix_2d
function matrix_A = assemble_matrix_2d(mesh, pde)
gauss_weights_tri = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];
dNdx = reference_basis_points_2d(mesh.basis_type, "dx");
dNdy = reference_basis_points_2d(mesh.basis_type, "dy");
ref_trans_tri = reference_basis_points_2d("P1");
matrix_A = sparse(mesh.Nb,mesh.Nb);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
    J_inv = inv(Jacobi);
    ref_points = [ref_trans_tri{1}(:), ref_trans_tri{2}(:), ref_trans_tri{3}(:)]*V;
    for i = 1:mesh.Nlb
        for j = 1:mesh.Nlb
            r = gauss_weights_tri * (pde.cor(ref_points(:,1),ref_points(:,2)) .* diag(([dNdx{i}; dNdy{i}]')*(J_inv*(J_inv'))*[dNdx{j}; dNdy{j}])) * mesh.absJ;
            matrix_A(mesh.elem(n,i), mesh.elem(n,j)) = matrix_A(mesh.elem(n,i), mesh.elem(n,j)) + r;
        end
    end
end
end

%% assemble_vector_2d
function vector_b = assemble_vector_2d(mesh, pde)
gauss_weights_tri = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];
psi = reference_basis_points_2d(mesh.basis_type);
ref_trans_tri = reference_basis_points_2d("P1");
vector_b = zeros(mesh.Nb,1);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    ref_points = [ref_trans_tri{1}(:), ref_trans_tri{2}(:), ref_trans_tri{3}(:)]*V;
    for j = 1:mesh.Nlb
        f = (pde.f(ref_points(:,1),ref_points(:,2))' .* psi{j});
        r = (gauss_weights_tri * f') .* mesh.absJ;
        vector_b(mesh.elem(n,j),1) = vector_b(mesh.elem(n,j),1) + r;
    end
end
end